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DIRECT  HUMERI  CAL  SIMULATIONS  OF  THE  PDF'S  OF  A  PASSIVE  SCALAR  IN  A 
FORCED  MIXING  LAYER 

P.  Givi  Flow  Research  Company,  21414-68th  Ave.  South,  Rent,  Washington  98032 

P.  A.  McMurtry  Department  of  Mechanical  Engineering,  University  of  Washington, 
Seattle,  Washington  98195 


Abstract  -  The  probability  density  functions  of  a  passive  scalar  quantity 
are  calculated  in  a  perturbed  mixing  layer  by  means  of  direct  numerical 
simulations.  The  results  indicate  that  the  two-dimensional  rollup  of  the 
unsteady  shear  layer,  and  the  pairing  process  in  particular,  contributes 
greatly  to  the  generation  of  the  predominant  peak  of  the  PDF's  within  the 
mixing  region. 

Key  Words  -  Direct  Numerical  Simulation,  Mixing  Layers,  Probability  Density 
Functions,  Coherent  Structures,  Entrainment. 


INTRODUCTION 

Probability  Density  Functions  (PDF's)  have  proven  very  useful  in  the 
theoretical  treatment  of  turbulent  reacting  flows  since  the  early  work  of 
Hawthorne  et  al.  (1949).  An  approach  based  on  the  solution  of  a  transport 
equation  governing  the  probability  density  function  of  the  scalar  quantities 
has  the  advantage  that  it  provides  a  complete  statistical  description  of  all 
the  scalars  (Pope,  1979).  Therefore,  the  effects  of  chemical  reactions  appear 
in  a  closed  form  eliminating  the  need  for  any  turbulence  modeling  associated 
with  the  scalar  fluctuations.  However,  models  are  needed  for  the  closure  of 
the  molecular  mixing  term  and  also  the  turbulent  convection  (O'Brien,  1981). 

In  most  of  the  previous  work  employing  the  PDF  approach,  the  effects  of 
molecular  mixing  have  usually  been  modeled  by  using  different  stochastic 
models  originating  from  the  same  "family"  of  the  coalesence-dispersion  models 
(Pope,  1982),  whereas  simple  gradient-diffusion  approximations  have  been 
employed  for  the  closure  of  the  turbulent  flux  of  the  PDF  (Givi  et  al.,  1984). 
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Among  these  previous  works,  Givi  et  al.  (1985)  used  a  Monte  Carlo 
numerical  routine  for  the  calculations  of  a  modeled  transport  equation 
governing  the  evolution  of  a  passive  scalar  PDF  In  a  nonreacting  two-stream 
turbulent  mixing  layer.  The  results  of  the  prediction  were  compared  with  the 
experimental  data  of  Masutani  and  Bowman  (1986),  which  were  obtained  under 
similar  hydrodynamlcal  conditions.  Good  agreement  between  predicted  results 
and  the  measured  data  was  obtained  for  the  first  two  moments  of  the  scalar 
quantity.  There  was  a  major  difference,  however,  between  the  calculated  and 
measured  profiles  of  the  PDF's.  The  experimental  results  indicated  that  the 
apparent  functional  form  of  the  PDF  changes  very  little  across  the  mixing 
layer  and  has  an  intermediate  peak  at  a  fixed  "preferred"  value  of  the 
concentration  (although  the  magnitude  of  this  peak  changes  as  the  mixing  layer 
is  traversed).  This  behavior  was  originally  documented  in  the  measurments  of 
Konrad  (1976)  and  Koochesfahani  (1984)  in  the  mixing  transition  and  post- 
mixing  transition  (Breidenthal,  1981)  regions  of  the  layer  and  Indicates  that 
a  given  mixed  fluid  concentration  has  the  same  probability  relative  to  other 
mixed  fluid  concentrations,  regardless  of  the  position  in  the  layer.  The 
predicted  results,  however,  indicate  that  the  location  of  the  PDF  peak,  with 
respect  to  the  concentration  coordinate,  changes  as  the  layer  is  traversed, 
meaning  that  the  PDF  of  the  mixed  fluid  concentration  varies  with  the  cross 
stream  direction  of  the  layer. 

The  major  reason  for  this  discrepancy,  as  suggested  by  Givi  et  al.  (1985) 
and  Masutani  and  Bowman  (1986),  is  due  to  the  shortcomings  associated  with  the 
gradient  diffusion  modeling  of  the  turbulent  flux  of  the  PDF  and  is  fairly 
Independent  of  the  modeling  of  the  molecular  mixing  term  (Kosaly  and  Givi, 
1987).  In  a  highly  intermittent  flow  such  as  a  mixing  layer,  regions  of 
turbulent  fluid  are  interrupted  by  the  presence  of  nonturbulent  surrounding 
fluid.  A  simple  gradient  diffusion  model  is  not  expected  to  accurately 
account  for  this  discontinuity. 

By  the  use  of  direct  numerical  simulations,  it  is  now  possible  to  simulate 
the  mixing  layer  directly  without  resorting  to  turbulence  models  (Riley  and 
Metcalfe,  1980).  Direct  numerical  simulation  refers  to  the  numerical  solution 
of  the  exact  transport  equations  of  turbulent  flows  by  means  of  very  accurate 
and  efficient  numerical  methods.  Transport  coefficients  are  chosen  to  assure 
that  all  relevant  flow  characteristics  are  accurately  resolved  so  that  no 
turbulence  modeling  is  required.  The  results  of  such  simulations  can  be  used 
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to  obtain  useful  information  on  the  statistics  of  the  variables  characterizing 
the  structure  of  the  flow,  which  in  turn  can  be  used  as  a  basis  for  turbulence 
modeling. 

In  this  communication,  the  results  of  direct  numerical  simulations  are 
used  to  construct  the  PDF  of  a  passive  scalar  quantltiy  in  a  perturbed  two- 
dimensional  mixing  layer.  The  profiles  of  the  PDF  constructed  in  this  manner 
are  used  to  address  the  shortcomings  associated  with  the  modeling  of  the 
turbulent  flux  of  the  PDF  in  the  transport  equation  governing  its  evolution. 

RESULTS 

A  pseudospectral  numerical  code  developed  by  McMurtry  et  al.  (1986)  was 
modified  to  simulate  a  two-dimensional  temporally  evolving  mixing  layer  under 
the  influence  of  harmonic  forcing.  The  numerical  resolution  was  upgraded  from 
the  previously  used  64  x  64  grid  to  256  x  256  equally  spaced  Fourier  modes. 
This  upgrade  was  required  for  better  statistical  analysis  of  the  data  used  for 
constructing  the  PDF's.  The  flow  is  assumed  periodic  in  the  streamwise 
direction  and  free  slip,  impermeable  boundary  conditions  are  employed  at  the 
transverse  boundaries.  Although,  the  laboratory  splitter  plate  mixing  layers 
evolve  spatially  downstream  and  the  numerical  simulations  evolve  temporally, 
Important  similarities  in  the  dynamics  of  these  two  flows  make  it  useful  to 
study  accurate  numerical  simulations  of  the  temporally  growing  layers.  By 
simple  Galilean  transformation,  a  flow  quantity  averaged  in  the  streamwise 
direction  can  be  related  to  the  time  average  of  the  same  quantity  at  a  fixed 
location  in  a  splitter  plate  configuration.  These  averaged  quantitities  are 
dependent  on  the  tranverse  coordinate  and  the  time.  Again,  the  inverse 
Galilean  transformation  relates  time  to  the  streamwise  location  in  a  splitter 
plate  configuration.  There  is  one  structure  within  the  periodic  domain  at 
each  time  step.  Statistical  analysis  are  perfromed  by  sampling  256  data 
points  in  the  streamwise  direction  at  each  transverse  location  and  at  each 
time  step.  The  presence  of  periodic  boundary  conditions  allows  us  to  use 
accurate  pseudospectral  numerical  methods;  these  methods  are  discussed  by 
McMurtry  (1987)  and  will  not  be  repeated  here. 

The  flow  field  is  initialized  with  a  hyperbolic  tangent  mean  streamwise 
velocity  profile  and  perturbations  corresponding  to  the  most  unstable  mode  of 
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this  profile  and  the  first  subharmonic  of  the  most  unstable  mode.  The 
properties  of  these  modes  have  been  evaluated  from  linear  stability  theory 
(Mlchalke,  1964).  The  fundamental  mode  in  the  mixing  layer  produces  a  single 
vortex  rollup.  When  the  subharmonic  is  added  in,  a  second  rollup,  or  pairing, 
can  occur. 

The  normalized  initial  value  of  the  conserved  scalar  concentration,  f, 
varies  from  0  in  the  bottom  stream  to  1  in  the  top  stream,  and  its  shape  is 
approximated  by  the  following  functional  form: 

y 

f(x,y,0)  -  -~j=  f  exp(-$/y  )2d* 

y0  /  T  J  u 
“OP 

The  flow  is  conveniently  characterized  by  two  nondimensional  parameters:  the 

Reynolds  number.  Re  ■  aOa/v,  based  on  the  mean  velocity  difference  across  the 

layer,  the  velocity  half  width,  and  the  kinematic  viscosity;  and  the  Peclet 

number,  Pe  »  Re  Sc,  where  Sc  is  the  molecular  Schmidt  number.  The  values  of 

these  two  parameters  were  set  equal  to  200,  so  that  the  scales  would  be 

accurately  resolved  on  the  256-256  grid  points  employed  in  the  numerical 

simulations.  The  value  of  y  was  chosen  so  that  the  initial  concentration 

o 

thickness  and  the  initial  velocity  thickness  were  identical.  The  time  depen¬ 
dent  transport  equations  governing  the  hydrodynamlcal  variables  (velocities 
and  pressure)  and  the  scalar  variable  (f)  are  solved  with  use  of  the  pseudo- 
spectral  code,  the  results  of  which  are  discussed  next. 

The  contour  plots  of  the  conserved  scalar  variable  are  shown  for  the 
purpose  of  flow  visualization.  In  Fig.  1,  we  present  the  time  development  of 
the  contours  of  the  variable  f  at  four  different  computational  times. 

Initially,  the  perturbation  associated  with  the  fundamental  mode  grows  until  a 
time  of  t*  (t*«tAu/s)  equal  to  12,  when  the  first  rollup  occurs.  Proceeding 
in  time  results  in  the  diffusion  of  the  core  of  the  vortex  and  the  growth  of 
the  subharmonic  mode,  which  expresses  Itself  in  the  form  of  a  second  rollup  and 
the  pairing  of  two  neighboring  vortices.  This  second  rollup  is  completed  by  a 
time  of  t*«36.  Proceeding  further  in  time  results  in  the  diffusion  of  the 
vortex  core  with  no  additional  rollup. 

The  profiles  of  the  PDF’s  of  the  variable  f  obtained  by  statistical 
analysis  of  the  instantaneous  values  of  f  are  shown  in  Fig.  2.  In  this  figure, 
the  PDF  has  been  plotted  as  a  function  of  the  instantaneous  concentration  f. 
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and  the  cross-stream  direction  of  the  shear  layer.  These  PDF’s  represent 
conditions  at  the  initial  stages  of  the  growth  (t*-3)  and  the  final  stage 
after  the  pairing  process  is  completed  (t*-36). 

A  comparison  between  the  two  parts  of  Fig.  2  reveals  the  effects  of  vortex 
dynamics  on  the  structure  of  the  PDF's.  Fig.  2a  shows  that,  at  the  early 
stages  of  the  development,  before  any  rollup  or  pairing  has  occurred,  the 
PDF's  can  be  simply  characterised  by  two  delta  functions  that  are  located  at 
f«0  and  f-1,  with  only  a  negligible  amount  of  mixed  fluid  at  the  mixing  core 
of  the  layer.  This  indicates  that,  at  any  location  in  the  cross-stream 
direction,  the  fluid  originates  from  either  the  top  stream  or  the  bottom  stream 
without  any  significant  amount  of  mixing  in  the  core.  Fig.  2b  shows  that  after 
the  occurance  of  the  rollup  and  the  completion  of  pairing,  a  third  spike 
appears  in  the  profiles  of  the  PDF  at  a  region  in  the  neighborhood  of  a  concen¬ 
tration  value  of  0.5.  The  presence  of  this  third  peak  suggests  the  existence 
of  a  preferred  mixed  fluid  concentration  equal  to  0.5.  The  combined  effects 
of  vortex  rollup  and  diffusion  are  to  engulf  fluids  from  the  two  streams  and 
mix  them  in  the  cores  of  the  vortices.  As  the  layer  is  traversed,  the 
preferred  value  of  this  concentration  does  not  seem  to  change  considerably  and 
remains  in  close  proximity  to  f«0.5. 

The  trlmodal  shape  of  the  PDF  is  consistent  with  that  observed  experimen¬ 
tally  (Roochesfahanl,  1984).  However,  the  experimental  measurements  were 
performed  in  the  transition  and  post-transition  regions  of  the  mixing  layer, 
whereas  the  numerical  data  presented  here  resulted  strictly  from  a  two- 
dimensional  simulation.  This  Indicates  that  the  two-dimensional  rollup  of  the 
unsteady  shear  layer  can  be  considered  as  one  possible  mechanism  by  which  the 
third  peak  is  generated  in  the  PDF  profiles. 

It  should  be  mentioned  that  the  presently  calculated  third  peak  of  the  PDF 
at  the  preferred  mixed  fluid  concentration  of  f-0.5  is  less  pronounced  than 
that  observed  experimentally  by  Masutanl  and  Bowman  (1986).  This  may  be  due 
to  the  fact  that  in  the  present  calculations,  the  effects  of  random  turbulent 
motion,  which  can  further  enhance  mixing,  are  not  taken  into  account.  In  order 
to  consider  the  contributions  of  turbulent  motion  into  the  generation  of  the 
third  peak  in  the  PDF  profile,  three-dimensional  simulations  are  required. 
Furthermore,  the  asymmetry  of  the  mixing  mechanisms,  which  has  been  both 
experimentally  (Roochesfahanl  et  al.,  1985;  Roochesfahanl  and  Dlmotakls,  1986; 
Mungal  and  Dlmotakls,  1984)  and  numerically  (Givi  and  Jou,  1987)  observed  for 
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spatially  evolving  flows  cannot  be  represented  In  the  temporal  simulations 
presented  here.  The  asymmetric  mixing  mechanism  results  in  a  preferred  concen¬ 
tration  value  that  is  closer  to  that  of  the  high  speed  stream  and  may  depend 
on  the  ratio  of  the  free  stream  velocities  (Dimotakis,  1986).  In  the  present 
temporal  simulations,  the  ratio  of  the  magnitude  of  the  free  stream  velocities 
is  equal  to  unity  and  the  structure  of  the  flow  is  symmetric  with  respect  to 
the  streamwise  coordinates  (as  shown  by  the  symmetry  of  the  PDF's  with  respect 
to  the  location  (f ,y)«(0.5,0)  in  Fig.  2).  Therefore,  the  numerical  simulations 
cannot  predict  any  other  than  the  arithmatic  average  of  the  concentration 
values  of  the  two  streams  (i.e.,  0.5).  Nevertheless,  the  results  indicate  that 
the  rollup  of  the  unsteady  shear  layer  contributes  greatly  to  the  generation  of 
the  predominant  peak  of  the  PDF's.  The  exact  role  of  the  small  scale  turbu¬ 
lence  motion  on  the  enhancement  of  such  generation  requires  full  three- 
dimensional  simulation  and  is  the  subject  of  our  future  research. 

Finally,  the  reason  that  the  methods  based  on  simple  gradient  diffusion 
modeling  of  the  turbulent  flux  of  the  PDF,  such  as  that  employed  by  Givi 
et  al.  (1985),  cannot  predict  the  trlmodal  shape  of  the  PDF  obtained  in  the 
present  simulation  is  due  to  the  fact  that  such  methods  do  not  consider  the 
influence  of  lntermlttency  caused  by  the  large  coherent  structures  in  the 
formulation.  The  results  of  the  direct  numerical  simulations  reported  here 
indicate  the  importance  of  such  structures  in  the  mixing  region  of  the  shear 
layer  and  suggest  the  need  for  better  turbulence  models  in  order  to  accurately 
predict  the  mechanisms  of  mixing  and  entrainment  in  such  flows. 
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FIGURE  CAPTIONS 


Figure  1 


Figure  2 


5642R 


:  Plot  of  the  Conserved  Scalar  Variable  (f)  Contours  at  Four 
Different  Computational  Times.  Contour  Minimum  is  0,  Contour 
Maximum  is  1,  Contour  Interval  is  0.1.  (a)  t*«3,  (b)  t*«12,  (c) 

t*-24,  (d)  t*-36. 

PDF’ 8  of  the  Conserved  Scalar  Variable  (f)  at  Points  Across  the 
Mixing  Layer  Width,  (a)  t*-3,  (b)  t*«36. 
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